Formation of longitudinal structures in granular flows 
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In the framework of the theory of partially fluidized granular flows we study the formation of 
longitudinal structures observed experimentally by Forterre and Pouliquen in a flow down a rough 
inclined plane. We show that the formation of longitudinal structures is related to the positive feed- 
back between the fluidization rate and the lateral stress (side pressure), which leads to a convective 
instability. Our theory explains main experimental features, such as appearance and amplification 
of the structure at some distance from the outlet and non-stationary behavior of the structures. 

PACS: 45.70.-n, 45.70.Ht, 45.70.Qj, 83.70.Fn 
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The dynamics of granular media has been an active 
area of research for physicists M and engineers @. One of 
the most interesting phenomena pertinent to many gran- 
ular systems is the transition from a static equilibrium to 
a granular flow. There has been a number of experimen- 
tal studies of flows in large sandpiles pH as well as in 
thin layers of grains on inclined surfaces |5|- 8| . Recent ex- 
periments with granular flows on a rough inclined plane 
revealed a new striking phenomenon: formation of non- 
stationary longitudinal structures at some distance down- 
stream from an outlet || . Authors of ^ proposed an ex- 
planation of the instability mechanism based on analogy 
with thermal convection in fluids: a rapid shear granu- 
lar flow leads to the increase of the granular temperature 
near the rough bottom. Because of the intrinsic dissi- 
pative nature of collisions between particles, the gran- 
ular temperature should decay away from the bottom, 
creating necessary conditions for the "thermal granular 
convection" (cf. |H]|ll]] ) . Using hydrodynamic equations 
obtained from the kinetic theory for dilute granular gases, 
this instability was studied analytically and numerically 
H . Although shear flow activated thermal granular con- 
vection could be a useful concept for the interpretation of 
experimental results, the theory || did not address two 
important observations: (1) experiments as well as sim- 
ulations [jl2l show that fluctuations of velocity are more 
significant near the free surface of the granular flow, and 
so the granular temperature may in fact be higher at the 
top rather than at the bottom, (2) longitudinal structures 
appear at some distance downstream from the outlet and 
exhibit complex spatio-temporal dynamics. 

In this Letter we demonstrate that the observed trans- 
verse instability can be described within the framework 
of the continuum theory of partially fluidized flows pro- 
posed by us in Refs. p~3)JT^ ] without invoking the con- 
cept of granular temperature more appropriate for dilute 
gas- like granular flows. In this theory, we introduce an 
equation for the order parameter which characterizes the 
phase state of the granular matter. In a certain range of 
parameters, the shear flow described by these equations, 



is unstable with respect to transverse perturbations. We 
show that this instability is of convective nature, i.e. 
small perturbations grow downstream while remaining 
small in the laboratory frame. Thus, the "rolls" appear 
at some finite distance from the outlet. This distance 
depends on the noise level and flow conditions. Since 
the pattern structure is determined by random fluctu- 
ations near the outlet, the resulting pattern is always 
non-stationary in the laboratory frame, very similar to 
that observed experimentally. 

Model. The starting point of our theory is the momen- 
tum conservation equation 



p Dvi/Dt 



Pogi, j = 1,2,3. 



(1) 



where Vi are the components of velocity, po = const is the 
density of material (in the following we set po = 1), g is 
acceleration of gravity, and D/Dt — d t + Vid Xi denotes 
the material derivative. Since the relative density fluctu- 
ations are small, the velocity obeys the incompressibility 
condition V • v = 0. The stress tensor is represented as a 
sum of the hydrodynamic part proportional to the flow 
strain rate and the strain-independent part, 



V 



dvi dvj 



(2) 



where r\ has the meaning of the viscosity coefficient. At 
this point we introduce the order parameter (OP) p which 
describes the "phase state" of the granular matter: it 
varies from in the "liquid" phase to 1 in the "solid" 
phase. We interpret the OP as the relative local den- 
sity of static contacts among the grains. In the solid 
state (p — 1) the strain- independent part should coin- 
cide with the "true" static stress tensor a% for the im- 
mobile grain configuration in the same geometry, and in 
the completely fluidized state (p = 0) we should have 
a~?j = —TISij (II is the hydrodynamic pressure). Accord- 
ing to our assumption, the off-diagonal elements of the 
strain-independent part of the stress tensor obey the con- 
ditions crfj — pa^j for i ^ j. In the solid state the normal 
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stresses a\ do not in general coincide. For the weakly- 
fiuidized state (p — ► 1) the normal stresses are close to 
the static values, however some dependence on the order 
parameter p (i.e. degree of fluidization) may appear. In 
the first order in 1— p one writes for the diagonal elements 
of the stress tensor 



al+a l {\-p) + 0{(l-pf) 



(3) 



where a, characterizes the response of the normal stresses 
on fluidization. Since fluidization is accompanied by the 
decrease in the number of static contacts among granules 
(i.e. dilution), the normal stresses should decrease with 
fluidization, i.e. oti > 0. It also agrees with the observa- 
tion in Ref. || that the crests of the surface deformations 
correspond to a more dilute granular state. 

According to Ref. Jl3|,0, the equation for the order 
parameter p is taken in the form 



p + vWp = V 2 p + p(l-p)(p-5), 



(4) 



where 5 = (</) — <fro)/((j)i — 4>o) is the control parameter, 
d> = max | <Jm n /Un n | is the maximum ratio of shear to 
normal stresses in the bulk, and tan -1 ifii^ are static and 
dynamic repose angles characterizing the granular mate- 
rial. Parameter cp which enters the Mohr-Coloumb yield 
condition Q , in the context of our theory is equivalent to 
a melting temperature in the theory of phase transition. 

Following the analysis of Ref. |l3| , we consider a layer 
of dry cohesionless grains on an inclined rough surface 
(see Fig.|l|). However, now we assume that the layer 
thickness h can vary in both x and y directions. The 
momentum conservation equation (|l|) in the stationary 
regime yields the force balance conditions 



<Jxz.x + Vyz,y + <J Z z,z = Q COS if, 
0~xx,x + °~xy,y + &xz,z = -gSmip 
0~xy,x ~t~ ®~yy.y ~t~ ®~yz.z 0? 



(5) 
(6) 
(7) 



where the subscripts after commas mean partial deriva- 
tives, and z = corresponds to the bottom of the layer. 
We assume that variations of the layer thickness along the 
y direction are small. Accordingly, there will be a small 
component of velocity along y direction and correspond- 
ing stress components a yz 

problem considered in Ref. |L3j these terms were irrele- 
vant. We show later that the transverse flux is crucial 
for the explanation of longitudinal structures. 
In the first order in h x ,h y Eqs. yield 




FIG. 1. Schematic representation of a chute geometry. 
2-axis is normal to the chute bottom, dashed line is paral- 
lel to the direction of gravity. 

For the chute geometry Fig. [j] parameter S in Eq. (jij) 
can be specified. For h = const there is a simple relation 
between shear and normal stresses, o~ xz — —t&n(po~ zz - 
The most "unstable" yield direction is parallel to the 
inclined plane, i.e. 4> — Wxz/czz\ an d 5 = So = (tan tp — 
</>o)/(</>i — ^o)- If h is a slowly varying function of x and 
y, one obtains for the parameter 8: 



5 = 8 Q -pK + 0{hi,hl) 



(9) 



where (3 — l/(0i — ^o) (see pSJTjl for more details). 

Thin layer solutions of Eq. (14) for the chute geometry 
Fig. [I] are subject to the following boundary conditions 
(BC): no-flux condition p z — at the free surface z = h, 
and no-slip condition p = 1 at the bottom of the chute 
z = (a granular medium is assumed to be in a solid 
phase near the rough surface). The velocity profiles cor- 
responding to a stationary profile of p{z), can be easily 
found from Eq. (|),(|), 



<1 



dv x 
dz 

OVy 

dz 



-g sin tp(z — h) 

CtydyPClZ 



-gs\rnp(l - p)(z - h), 
(10) 



and (Jyy In the avalanche In thc deriv ation of Eq. (|1C|) we used = for a flat 

layer. The components of the flux J = J Q v(z)dz are 



Jx 



g simp 



dz / (1 - p(z'))(z' - h)dz' 
V Jo Jo 



<7zz ~ g cos (p(z - h) , o xz w -gsia<p(z - h), 



J a 



a 



/ dz 
V Jo 



dz" / dypdz' 



(11) 



dza 



yy,y 



(8) 



Here we used the conditions a xx — const, o~ xy = const 
and <j zz — a xz = a yz = on free surface z = h. 



For thin weakly- fluidized layers (see [|13 14|), we can 
look for solution of Eq. (0) in the form 



1 — A sin 



m — z 
\2h J 



h.o.t., 



(12) 
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where A C 1 is a slowly varying function of t, x, and 
y. Substituting ansatz (O) into Eq. (^) and applying 
orthogonality conditions ]13| , P^ |), we obtain 

A t = AA + X7\A + 8(2 ~ 6) A 2 - ^A 3 - fih 2 Ad x A (13) 
3tt 4 

where Vi = 8* + d%, A = 8 - 1 - p, = (3tt 2 - 
16).gsiniy9/37r 3 77 = 0.146g sin ip/rj. The mass conservation 
yields 

dh 

-jfi = -V • J = ~d X J X - dyJy (14) 

where the flux components J x<y calculated from Eq. jll| ) 
in the thin layer approximation, are 

J x = fih 3 A , J y = -aih 2 Ah y + a 2 h 3 A y (15) 

/1 = 2(tt 2 — 8)gsm.tp/r)Tr 3 w 0.12g sin (p/rj, a\ = o- y (^ 3 + 
8tt - 48)/2tt 3 77 w 0.13a y /r/, a 2 = 8a y /ir 3 ri w 0.26a y /rj. 
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FIG. 2. Phase diagram in S — Jo plane. 

Our model has two control parameters, the downhill 
mass flux Jo at the outlet x = 0, and 5. For large 
enough Jo and S, Eqs. (|l^),(|l^) possess a steady-state 
solution A = Ao,h = ho (see Fig.|J). The stability of 
this solution can be examined by substituting the ansatz 
{A, h} — {Aq, ho} + {£, £} exp(Ai + ik x x + ik y y), where 
{C)£} is a small perturbation. For long- wave perturba- 
tions (k x ,k y «C 1) in the leading order one obtains for 
the growth rate A 

X = ik xfJ ,hlA Q + Bk 2 y + ... (16) 

where B = 7r 2 a 2 [16(2- ( 5)/37r-3 J 4o] _:L -ai/ioAo. As seen 
from Eq. (JTsj) , pure longitudinal perturbations k y = 
are neutrally stable, but for B > the perturbations 
with transverse component of the wave number can be 
unstable. For pure transverse perturbations (k x — 0) the 
instability is aperiodic (A is real) , as in thermo-convection 
in ordinary fluids. The physical meaning of the instabil- 
ity is the following: decrease of p (increase of fluidization) 
corresponds to the decrease of the lateral pressure {—cr yy ) 
which in turns lead to the increase of the height h and 



further fluidization. The phase diagram of the instabil- 
ity on S — Jo plane is shown in Fig. ^. The transverse 
instability exists in the band restricted from above by 
the condition B = and from below by the existence of 
steady-state solution |Q . 

Complexity of A and asymmetry of the problem in re- 
direction signal the possibility of the convective nature of 
the instability: the perturbations may grow in the mov- 
ing frame and decay in the laboratory frame (see, e.g. 
|l5],[l6|] ) . As a test for convective instability one has to 
examine the value of A for the saddle point of the func- 
tion \(k x ,k y ) in the complex k x plane, i.e. d\/dk x = 0. 
In fact, it is easy to see that if the instability is weak (for 
a y <C 1) and the downhill flux J is finite, the instability 
has to be convective. In the presence of persistent fluc- 
tuations, e.g. at the flow outlet, the perturbation will be 
spatially amplified down the flow. 

Let us estimate the spatial growth rate. Consider the 
unstable perturbations caused by a small stationary y- 
periodic force with wavenumber fco localized at x = 0. 
Calculations show that the linearized solution £ will be 
described by the integral 

r°° dk 
C~ y_ exp[ l kr]^y(l-exp[A(k)t]) (17) 

The non-stationary part ~ exp[iA(k)r] decays in time 
due to the convective character of the instability, and 
the integral yields 

C ~ exp[ifc y + ik*x] (18) 

where k* is found from X(k x — k*, kg) — 0. It is easy to 
see from Eq.(|i~6|) that, for small fc°, k* « ~iBk 2 /fih 2 A . 
Since k* is imaginary, perturbations from a stationary 
source at x — grow exponentially downhill. For random 
perturbations introduced at the outlet the typical scale 
of the pattern will be determined by the most unstable 
wavenumber in the transverse direction. However, the 
pattern will remain non-stationary due to intrinsically 
random nature of the noise. 

We studied Eqs. (r^),(E3) numerically. The simula- 
tions were performed in a large system, more than 2000 
dimensionless units in x-direction (downhill), and 200 
units in y-direction. Fixed flux boundary conditions were 
imposed at the outlet at x — 0. Selected results are 
presented in Figs. ||^. Figure || illustrates spatial am- 
plification of perturbations downstream and formation 
of longitudinal structures. As it follows from our simula- 
tions, far away from the outlet the structures remain non- 
stationary and exhibit spatio-temporal dynamics which 
are very similar to observed in Ref. ||. The profiles of 
h and A vs y are shown in Fig. ||. In agreement with 
experiment, crests in h corresponds to crests in A, i.e. to 
more fluidized regions of flow. 
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FIG. 3. Grey-scale snapshots of h(x, y) (white corresponds 
to larger h), for S = 1.5, n = 0.025, (3 = 3.14, J = 1.75, 
a y — 0.04. Equilibrium value of the layer thickness ho ~ 4.45. 
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FIG. 4. Profiles of h(y) and A(y) at x = 2000. 

We demonstrated that formation of longitudinal struc- 
tures in granular flow down a rough inclined plane is the 
result of noise amplification and saturation due to the 
convective instability. The mechanism of the instability 
is related to the dependence of the local pressure on the 
fluidization rate controlled by the order parameter. We 
conjectured a simple linear relation between these quan- 
tities. It would be of interest to verify this relation in 
physical experiment and molecular dynamic simulations. 
The convective character of the instability suggests that 
the flow can be controlled by small systematic perturba- 
tions introduced at the outlet. 

In recent paper pcfl , a novel fingering instability in a 
thin granular layer inside a horizontal rotating cylinder 
was reported. We believe that the nature of this insta- 



bility fundamentally is the same as described here. The 
only important difference is that due to backflow the per- 
turbation can propagate upstream, and thereby the con- 
vective instability becomes "global". This may explain 
why this fingering instability occurs in a "short" system 
as compared with inclined layer experiment B . 
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